In this lab you will learn



library(tidyverse)
library(gridExtra) ## for plotting function grid.arrange()
library(deSolve) ## for predator-prey-resource ODE model
library(magrittr) ## for pipe symbol %<>%
library(knitr) ## for table-viewing function kable()
theme_set(theme_bw())
theme_update(
  panel.grid = element_blank(),
  aspect.ratio = 1
)



Parasite-host and pathogen-host interactions are among the most common and most important ecological interactions in nature. Virtually every organism on Earth has its pathogens – even some pathogens have their own pathogens! – and in many cases it is infectious disease, not predation or resource limitation, which keeps a population in check.

The ecology of infectious disease is not very different from that of predator-prey or consumer-resource dynamics: the pathogen plays the part of the predator/consumer, and the host is the prey/resource. However, there are important differences.


Figure: Bacteriophages attached to a bacterial cell, by Dr Graham Beards. Blood smear of Plasmodium falciparum by Lukas.S. Electron micrograph of SARS-CoV-2, by NIAID.



The SIR model

Those differences and similarities allow us to study infectious disease dynamics using similar-but-different models as the ones we saw in the previous module. In particular, the Susceptible-Infected-Resistant (SIR) model is the mainstay of disease ecology. This model makes the following assumptions:

We can put the words above in mathematical terms:

\[ \begin{eqnarray} \frac{dS}{dt} &=& -\beta SI \\ \frac{dI}{dt} &=& \beta SI - \gamma I - \mu I \\ \frac{dR}{dt} &=& \gamma I \end{eqnarray} \] where the symbols in the model have the following meanings

Symbol Meaning
\(S\) percentage of susceptible individuals in the population
\(I\) percentage of infected individuals in the population
\(R\) percentage of recovered individuals in the population
\(\beta\) infection rate
\(\mu\) disease-caused mortality
\(\gamma\) recovery rate

We can immediately draw some conclusions from the equations above:

\[\frac{dI}{dt} = \beta I\left(S - \frac{\gamma + \mu}{\beta}\right)\]

The R code for running and plotting the SIR model follows below:

SIR_Model = 
  function(
    initial_S,
    initial_I,
    initial_R,
    infection_rate,
    recovery_rate,
    mortality_rate,
    final_time,
    time_step
  ){
      SIR = 
        function(t, state, parameters){
           with(as.list(c(state, parameters)), {
             dSdt = -beta * S * I
             dIdt = beta * I * S - nu * I - mu * I
             dRdt = nu * I
             
             list(c(dSdt, dIdt, dRdt))
           })
          }
      
      times = seq(0, final_time, by = time_step)
      
      parameters = 
        c(
          beta = infection_rate,
          nu = recovery_rate,
          mu = mortality_rate
        )
      
      state = 
        c(
          S = initial_S,
          I = initial_I,
          R = initial_R
        )
      
      out = ode(y = state, times = times, func = SIR, parms = parameters)
      
      return(
        list(
          parameters = 
            list(
              infection_rate, 
              recovery_rate, 
              mortality_rate
            ),
          initial_conditions = list(initial_S, initial_I, initial_R),
          state = out
        )
      )
   }

Plot_SIR = function(model){
  as.data.frame(model$state) |>
    pivot_longer(-time, names_to = 'type', values_to = 'proportion') |>
    mutate(type = factor(type, levels = c('S', 'I', 'R'))) |>
    ggplot(aes(time, 100 * proportion)) +
    geom_line(aes(color = type), linewidth = 1) +
    facet_wrap(~type, scales = 'free') +
    labs(
      x = 'time (days)',
      y = 'percentage of the population'
    ) +
    theme(
      legend.position = 'none',
      strip.background = element_rect(fill = 'orange')
    )
}



\(R_0\) and the fate of the disease

Let’s first examine the simplest possible scenario, where there is no disease-induced mortality \(\mu=0\) and the recovery rate \(\gamma\) is zero or much lower than the infection rate \(\beta\). This corresponds to a case where the host has little to no immunity to the disease, possibly because it is a newly evolved pathogen that the host had previously never encountered. Let’s set \(\beta = 2\) and suppose 0.1% of the host population is initially infected.

What do you predict will happen in terms of the proportion of the population that is infected?

Plot_SIR(
  model = 
    SIR_Model(
      initial_S = .999,
      initial_I = .001,
      initial_R = 0,
      infection_rate = 2,
      recovery_rate = 0,
      mortality_rate = 0,
      final_time = 30,
      time_step = .1
    )
) +
  ggtitle('Scenario 1')

Perhaps unsurprisingly, the disease sweeps through the entire population.


Q1. What happens if we decrease the infection rate? Will the entire population still become infected? If so, what is different about that scenario?


Now let’s add a recovery rate, \(\gamma = 1\):

Plot_SIR(
  model = 
    SIR_Model(
      initial_S = .999,
      initial_I = .001,
      initial_R = 0,
      infection_rate = 2,
      recovery_rate = 1,
      mortality_rate = 0,
      final_time = 30,
      time_step = .1
    )
) +
  ggtitle('Scenario 2')

The disease still spreads initially, but the spread eventually slows down and halts, and the disease eventually fades away. This may sound great, but look what happened to the susceptible curve. By Day 30, less than 25% of the host population never contracted the disease. In other words, over three quarters of the population got sick!

Does this scenario qualify as an epidemic?

According to the CDC,

Epidemic refers to an increase, often sudden, in the number of cases of a disease above what is normally expected in that population in that area. Outbreak carries the same definition of epidemic, but is often used for a more limited geographic area.

In other words, whether this is an epidemic or not depends on what is normally the case with this disease. Again according to the CDC,

While some diseases are so rare in a given population that a single case warrants an epidemiologic investigation (e.g., rabies, plague, polio), other diseases occur more commonly so that only deviations from the norm warrant investigation.


Q2. Let’s assume there is no disease-caused mortality, i.e. \(\mu=0\), and consider two scenarios: a) high transmissivity and low recovery, i.e. infection rate \(\beta = 10\) and recovery rate \(\gamma = 0.1\); b) low transmissivity and high recovery, \(\beta = 0.1\) and \(\gamma = 10\). What proportion of the population eventually becomes infected in each case? (Hint: look at the recovered class.)


According to the SIR model, if you are infectious, the rate at which you infect another person per unit time is \(\beta\), whereas the rate at which you leave the infectious class (via either recovery or death) is \(\gamma + \mu\). Therefore, each currently infectious individual contributes a new individual to the infectious class (a new “birth” for the disease) at rate \(\beta\) and removes herself from the class (a “death” to the disease) at rate \(\gamma + \mu\). If the per capita rate of removals surpasses that of additions, the disease cannot grow. This suggests that whether or not the disease is capable of spreading depends on whether the ratio \(R_0 \stackrel{\text{def}}{=} \frac{\beta}{\gamma + \mu}\) is greater or less than 1. Indeed, epidemiologists recognize this important quantity as the basic reproductive number of the disease.


Q3. Calculate \(R_0 \equiv \frac{\beta}{\gamma + \mu}\) in each of the zero-mortality SIR scenarios above: a) No immunity \(\beta=2\), \(\gamma=0\); b) Some immunity \(\beta=2\), \(\gamma=1\); c) High transmissivity \(\beta=10\), \(\gamma=0.1\); d) Low transmissivity \(\beta=0.1\), \(\gamma=10\). Assuming each of these scenarios corresponds to a different pathogen, put those diseases in order of ability to spread. Which of them can become epidemic in a naive population?


Consider the point in time where the number of infections stops growing and starts falling. By definition, this is the point where the derivative \(\frac{dI}{dt}\) is zero. Looking back at the model equations, we conclude that on this day, the proportion of the population that is susceptible to the disease reaches a critical value

\[ S_c = \frac{\gamma + \mu}{\beta} \equiv \frac{1}{R_0} \]

In words, on the day when the infected class stops growing, the proportion of the population who has not yet contracted the disease will reach a critical value \(S_c=\frac{1}{R_0}\). This provides us with a way to empirically estimate the \(R_0\) from data on disease cases.


Figure: Herd immunity diagram by Narayana Health


Herd immunity

Look again at the equation for \(\frac{dI}{dt}\) and ask yourself what happens if \(S(t) < S_c\) in terms of the ability of the disease to spread any further.

We can see that \(S_c\) works as an important threshold: When the proportion \(S\) of susceptible individuals in the host population is lower than that threshold, i.e. when \(S < S_c\), the disease dies out.

Note: Notice that if \(R_0 < 1\), then the threshold for \(S\) is \(S_c > 1\). Since the maximum value \(S\) can take on is 1, this means that the condition \(S < S_c\) is automatically satisfied in this case and the disease can never spread. This is consistent with our interpretation of \(R_0\) as the basic reproductive number of the disease.

Consider everyone in the host population who are not currently in the susceptible category, \(1 - S\). Let’s call that proportion \(p\). We can rewrite our threshold \(S_c\) as a threshold for \(p\):

\[ p_c = 1 - S_c = 1-\frac{1}{R_0} \] How should we interpret \(p_c\)? Notice that if the proportion of the population who cannot catch the disease is at least as high as \(p_c\), the disease is starved of susceptible hosts and dies out. Accordingly, \(p_c\) is called the herd immunity threshold, or HIT.

Note: Note that \(p_c\) is interpreted as a proportion, and therefore must be contained within the range \([0, 1]\). Thus, the definition above is nonsensical if \(R_0<1\), and should be simply replaced with 0 in that case. So we should really be saying that

\[\begin{equation} \textrm{HIT} = \begin{cases} 1-\frac{1}{R_0} & \text{if}\ R_0 \geq 1 \\ 0 & \text{otherwise} \end{cases} \end{equation}\]


Q4. Calculate the herd immunity threshold in each scenario in Q3. Give the answer in percentages of the host population.


Who is in the category \(p\)? Everyone who cannot catch the disease, either because they already have it, or because they’ve had it and became immune, or … vaccines.

The core concept of herd immunity is that a disease can be eradicated even if there are still susceptible individuals in the population, i.e. even if only a percentage of the population is immunized.


Note

To answer Q5 and Q6 below you need to open the I herd you! app on your browser by clicking on this link.


We can get a visual sense of what herd immunity looks like by navigating to the I herd you! app . This represents a simulation of an SIS model (susceptible-infected-susceptible), which differs from an SIR model in that individuals who recover from the infection can be reinfected. The page lets you simulate multiple scenarios of contact between individuals: mixing population, static network, dynamics network, and lattice model. Let’s focus on the static network option, as the disease dynamics plays out faster in that scenario.

Upon clicking the Play button, you will see the spread of the disease across the population. Notice how the disease becomes endemic in the population, meaning there are always some individuals who are currently infected.


Q5. How do you explain the fact that the disease always eventually dies out in the SIR model but not in the SIS model?


Now let’s vaccinate some people. Vaccinated individuals become temporarily immune to the disease. Slide the Vaccine Uptake button to about half of its length, and watch what happens within.


Q6. What eventually happens? What does this show in terms of the need to vaccinate everyone for eradication of the disease?

Feel free to explore the other population scenarios. This app is called Explorables after all.


Q7. Would the flu be better described by a SIR model or an SIS model?


Figure: Vaccines affect disease epidemiology by removing susceptible individuals from the host population. Illustration from Modern Healthcare


Vaccines

As you just saw in answering the questions above, vaccines do double duty in fighting an epidemic. They keep you from getting sick, and they stop you from transmitting it to others.

Suppose at the start of an outbreak, an efficient public health campaign manages to get a proportion \(\varphi\) of the population vaccinated. Assuming the vaccine is 100% efficient, this brings the proportion of the initial population that can be infected down to \((1-\varphi)\). What effect will that have on the epidemiology of the disease?

Let’s think back on the components of the infection rate \(\beta\): this is the rate at which an infected individual infects healthy individuals. We can write it as the product of the contact rate \(c\) between people and the probability that the encounter leads to transmission. However, now that only a proportion \(1-\varphi\) of healthy individuals can contract the disease, we must write the new contact rate as \(c(1-\varphi)\), reflecting the lower rate at which infected individuals will meet susceptible individuals. In effect, we have a new infection rate \(\beta' = (1-\varphi)\beta\). What impact does that have on the prospects of an epidemic?

Essentially, we will have a new reproductive number, \(R_e = (1-\varphi)R_0\), called the effective reproductive number of the disease. Now that some people are immune, this becomes the important number to compare to 1. In other words, whether the disease spreads or dies out depends on whether \(R_e\) is greater or less than 1.


Q8. Consider a disease with epidemic capability – that is, its basic reproductive number \(R_0\) is greater than 1 – which has been just introduced to a population. We wish to vanquish this disease with a vaccine, which we will assume to be 100% effective. What is the critical proportion of the population that must be vaccinated to stop the disease?


Notice that by artificially reducing the susceptible class and enlarging the resistant class, successful vaccines essentially bring the host population forward in time to the herd immunity stage, without the human cost of widespread infections.

Of course, vaccines are rarely 100% effective, meaning they may not bring the risk of acquiring a disease down to zero.


Q9. In class we discussed vaccine efficacy. a) Explain in your own words what vaccine efficacy is and how it is estimated using clinical trials. b) Consider a hypothetical trial for a new vaccine against a certain disease. If 1% of volunteers who received a placebo contracted the disease, while only 0.05% of those who received the vaccine were infected, what is the efficacy of the vaccine?


According to the New York Times, “In 2020, the F.D.A. set a goal for coronavirus vaccine trials. Each manufacturer would need to demonstrate that a vaccine had an efficacy of at least 50 percent. The confidence interval would have to reach down no lower than 30 percent. A vaccine that met that standard would offer the kind of protection found in flu vaccines — and would therefore save many lives.”


Q10. Suppose a vaccine manufacturer seeking approval from the F.D.A. reports trial results as follows: Out of 20,000 participants, 10,100 were vaccinated and 9,900 received a placebo. In the vaccinated group, 10 individuals tested positive for covid within 28 days. In the placebo group, 20 individuals tested positive in the same period. a) Calculate the estimated efficacy of the vaccine. b) Use the code below to estimate the 95% confidence interval around that point estimate. Would this vaccine meet the F.D.A. approval standards?

CI95 = function(CV, NV, CP, NP){
  RR = (CV / NV) / (CP / NP)
  varlogRR = 1 / CV - 1 / NV + 1 / CP - 1 / NV
  CI_lower = 100 * (1 - RR * exp(+1.96 * sqrt(varlogRR)))
  CI_upper = 100 * (1 - RR * exp(-1.96 * sqrt(varlogRR)))
  return(round(c(CI_lower, CI_upper)))
}


Note

To answer Q11 through Q13 below you need to open the Facebooked Flu Shots app on your browser by clicking on this link.


One last thing to know about vaccines is that some mass vaccination strategies are more efficient than others. Let’s see how that works using the Facebooked Flu Shots app. This app will show you a network of people, which you can think of as friends who are in regular contact with each other, and therefore can transmit a disease to each other. Not everyone is friends with everyone else and some people have more friends than others. Just like in real life. A network therefore serves as a conduit for a disease, which can go from person A to person Z via their intermediate friends.

The idea here is that vaccines disrupt those conduits by removing some people from the transmission network (while hopefully keeping friendships intact). So the question is, if we are going to vaccinate a certain percentage of all the people, can we do any better than vaccinating a random sample of the population?

Open the app, click on the Erdos-Renyi network option (which makes prettier-looking networks), leave the vaccination percentage at 38%, and click on Option A: random vaccination. You will notice that a bunch of people get removed from any networks – either because they are not friends with the people who remain connected by networks, or because they were vaccinated, in which case they become isolated (from the disease’s perspective). Notice how even though we vaccinated about a third of the population, many unvaccinated people are removed from the networks and thus become insulated from the disease. That’s herd immunity at work. However, many other people still remain vulnerable, and indeed the disease spreads through those networks.


Q11. Now try Option B, strongly connected, which targets people with lots of friends to be vaccinated. The same number of people are vaccinated, but this time we see very different results. Explain the difference.


Finally, try Option C. We sample the same percentage of the population at random as in Option A, but this time we choose a random friend of each person sampled to get vaccinated.


Q12. Read the explanation on the app, then explain in your own words why Option C works better than Option A.


Q13. Having seen how targeted vaccination strategies can be more efficient than fully random ones, discuss some reasons why we may be stuck with Option A from a public policy perspective, even though it is the least efficient. You are welcome to speculate, but good places to look for information on this include newspapers, the World Health Organization, the National Institute of Health and other governmental websites, and Wikipedia (make sure to check their references.).



Bonus Questions (Optional)

One of the most important skills for a student of the health sciences is the ability to interpret a positive result on a test that screens for a disease. Suppose a 62-year old female patient is tested for breast cancer, which affects about 350 in 100,000 women in that age category in the US. The test is a recently developed blood-based test with sensitivity 80% and specificity 79%.

Q14. (1 extra point in the lab) If the patient tests positive, what is the probability she has breast cancer?

Q15. (1 extra point in the lab) If the patient tests negative, what is the probability she has breast cancer?

In order to answer these questions, we must understand how the different statistics above play into the probabilities we wish to know. And in order to do that, we must use Bayes’ theorem, which is a theorem about conditional probabilities. A theorem is a mathematical proposition which may not be self-evident but is demonstrably true. A conditional probability is the probability of something we wish to know—our hypothesis—given something that we do know—the evidence. In our case, we want to know the probability that the patient has cancer (the hypothesis) given that she tested positive (the evidence).

Before you are shown any evidence, you may have prior beliefs about the probability of the hypothesis. For example, in the case of breast cancer, we can say our patient has a 350 in 100,000 chance of having cancer even before she steps into the doctor’s office because that is the prevalence of the disease among women of her age.

Bayes’ theorem is all about updating our beliefs regarding our hypothesis once we discover some new evidence. In statistical lingo, Bayes’ theorem tells you how to update the prior probability of the hypothesis into a posterior probability of the hypothesis given the evidence. Specifically, it says that to arrive at the posterior (our new belief), you must multiply the prior (our old belief) by the likelihood of the evidence assuming the hypothesis is true, divided by the overall probability of the evidence. In our case, this corresponds to multiplying our prior belief about the hypothesis that the patient has cancer (350 in 100,000) by the likelihood that her test would come out positive if she did have cancer, and then dividing the result by the overall probability that the test would come out positive on a random person. This makes sense: the more likely it is that a person with cancer would get a positive test result, the more likely it is that a who tests positive has cancer. On the other hand, the more likely a test is to come out positive overall, the less confident we are that our patient has cancer just because she tested positive (after all, it could have been a false negative).

Let’s put the words above into a precise mathematical expression. If we use the notation \(P(A|B)\) to represent the conditional probability that proposition A is true given that we know proposition B is true, then Bayes’ theorem reads

\[ P(A|B) = \frac{P(B|A)}{P(B)}P(A) \] This expression may look opaque, but it says the exact same thing as the paragraph above. In particular, if we identify proposition A with the hypothesis (cancer) and propositon B with the evidence (positive test result), and use the names posterior, prior, and likelihood as explained above, we can rewrite the theorem as follows:

\[ \textrm{Posterior}(\textrm{cancer}|+) = \frac{\textrm{Likelihood}(+|\textrm{cancer})}{\textrm{Prob}(+)} \, \textrm{Prior}(\textrm{cancer}) \]

We’re almost ready to plug in the statistics from above about disease prevalence and the properties of the test. \(\textrm{Posterior(cancer|+)}\) is what we wish to know, \(\textrm{Prior(cancer)}\) is our prior belief about the probability our patient has cancer (350 in 100,000, or 0.35%), and \(\textrm{Likelihood(+|cancer)}\) is the sensitivity of the test (80%). But we still need to figure out the overall probability that the test would come out positive on a random person, \(\textrm{Prob(+)}\).

There are two ways for a person to test positive: either they do have cancer and the test correctly identifies it (a true positive), or they do not have cancer and the test falsely indicates that they do (a false positive). A true positive will occur with probability \(\textrm{Likelihood(+|cancer)} \cdot \textrm{Prior(cancer)}\); that is, we multiply the prior probability that a person has cancer by the likelihood that a person with cancer would test positive. Analogously, the second possibility (false positive) will occur with probability \(\textrm{Likelihood(+|no cancer)} \cdot \textrm{Prior(no cancer)}\); that is, we multiply the prior probability that a person does not have cancer (which is 100% minus the probability that the person does have it) by the likelihood that a person without cancer would test positive.

The last step is to recognize that \(\textrm{Likelihood(+|no cancer)} + \textrm{Likelihood}(-|\textrm{no cancer}) = 1\). In words, if we do not have cancer, we either will test positive or negative with total probability 100%. This is convenient because it allows us to write \(\textrm{Likelihood(+|no cancer)} = 1 - \textrm{Likelihood}(-|\textrm{no cancer})\). Note that \(\textrm{Likelihood}(-|\textrm{no cancer})\) is the specificity of the test—i.e. the probability that it returns negative when we do not have the disease.

Putting these things together gives us

\[ \textrm{Prob(+)} = \textrm{Likelihood(+|cancer)} \cdot \textrm{Prior(cancer)} \; + \; \left(1 - \textrm{Likelihood}(-|\textrm{no cancer})\right) \cdot \left(1 - \textrm{Prior(cancer)} \right) \]

You can now answer Q13 by plugging the right statistics in the right places, and after a little bit of thought, you should also be able to figure out how to apply Bayes’ theorem to answer Q14.

One lesson to draw from all this is that we can’t know the probability of having a disease upon testing positive for it if we don’t know the prevalence of the disease in the population unless there is 0% chance of a false positive. Testing positive on a test that correctly identifies the disease in 99% of those who have it emphatically does NOT mean that you have 99% chance of having the disease. The true probability is lower, and for rare diseases it can be much lower.

Beyond testing for diseases, Bayes’ theorem is critical in any situation where we need to know how to update our beliefs based on evidence. This includes policymaking, public health planning, and criminal trials. Imagine how easy it would be to misjudge the probability of guilt of a defendant in the face of incriminating evidence if the jury or judge have a poor grasp of the prior probability of guilt, the probability that the evidence may be faulty (for example, false identification by a mistaken witness), and how all this comes together in Bayes’ theorem! Bayesian thinking has enormous potential to avoid miscarriage of justice in courtrooms but unfortunately seems to have consistently met with resistance from the courts. My favorite quote from the blog post linked above is as follows:

There is a persistent attitude among some members of the legal profession that probability theory has no role in the courtroom. In one case in England, in fact, an appeals court denounced the use of Bayesian calculations, asserting that members of the jury should apply “their individual common sense and knowledge of the world” to the evidence presented.

Apart from the obvious idiocy of using common sense to resolve complex issues, the court’s call to apply “knowledge of the world” to the evidence is exactly what Bayesian math does. Bayesian reasoning provides guidance for applying prior knowledge properly in assessing new knowledge (or evidence) to reach a sound conclusion. Which is what the judicial system is supposed to do.


\[ \color{red}{***\; \textrm{This is the end of Part 1 of Lab 7} \;***} \]